import numpy as np
import matplotlib.pyplot as plt # plt 用于显示图片
import matplotlib.image as image # mpimg 用于读取图片
from skimage import morphology
from skimage.color import rgb2gray
from sklearn.model_selection import cross_val_score
from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor
from sklearn import preprocessing
from sklearn import metrics

from com.csu.xubobo.feature_extract import feature_extract1


def findlabel(label_image, x, y):
    r_label = np.copy(label_image[x - 2:x + 3, y - 2:y + 3, 0])
    b_label = np.copy(label_image[x - 2:x + 3, y - 2:y + 3, 2])
    if r_label.sum() > b_label.sum():
        return 1
    elif r_label.sum() < b_label.sum():
        return 0
    else:
        return 2
x_data=[]
y_data=[]
LOCATION=[]
#读入图片
for num in range(1,21):
    print(num)
    img = image.imread('/Users/xubobo/Desktop/毕业设计/代码/各种图片/DRIVE/test/images/'+str(num)+'_test.tif') # 读取和代码处于同一目录下的 lena.png
    #print(img.shape) #(584, 565, 3)
    plt.imshow(img);plt.axis('off');plt.show()
    #print(img.dtype)
    segment= image.imread('/Users/xubobo/Desktop/毕业设计/代码/各种图片/DRIVE/test/1st_manual/'+str(num)+'_manual1.gif',format='grayscale')    #血管分割结果
    segment = rgb2gray(segment)
    #print(segment.shape) #(584, 565)
    #plt.imshow(segment);plt.axis('off');plt.show()

    mask = image.imread('/Users/xubobo/Desktop/毕业设计/代码/各种图片/DRIVE/test/mask/'+str(num)+'_test_mask.gif')          #视野
    #print(mask.shape) #(584, 565)
    #plt.imshow(mask);plt.axis('off');plt.show()

    od = image.imread('/Users/xubobo/Desktop/毕业设计/代码/各种图片/DRIVE/GrowSegment1/'+str(num)+'_test_seg.png')          #视盘
    #print(od.shape) #(584, 565)
    #plt.imshow(od);plt.axis('off');plt.show()

    label_image = image.imread('/Users/xubobo/Desktop/毕业设计/代码/各种图片/DRIVE/DRIVE（动静脉分类金标准图像）/'+str(num)+'_test_ManualAV.tif')          #视盘
    #print(label_image.shape) #(584, 565,3)
    #plt.imshow(label_image);plt.axis('off');plt.show()
    r_label=np.copy(label_image[:,:,0])
    g_label=np.copy(label_image[:,:,1])
    b_label=np.copy(label_image[:,:,2])
    label_image=np.zeros(label_image.shape)
    label_image[:,:,0]=r_label-g_label
    label_image[:,:,2]=b_label-g_label
    label_image[label_image>0]=1
    #plt.imshow(label_image);plt.axis('off');plt.show()

    #预处理

    #提取感兴趣区域
    odcenter,oddis=morphology.medial_axis(od,return_distance=True)      #获取视盘中心和半径
    #plt.imshow(odcenter);plt.axis('off');plt.show()
    x_temp,y_temp=np.nonzero(odcenter)
    x_od=np.mean(x_temp)
    y_od=np.mean(y_temp)
    oddistance=np.max(oddis)
    #print(x_od)
    #print(y_od)
    #print(oddistance)
    #


    roi=np.zeros(od.shape)      #感兴趣区域
    ir,ic = np.indices(roi.shape)
    circle1 = (ic - y_od)**2 + (ir - x_od)**2 < (oddistance * 3)**2
    circle2 = (ic - y_od)**2 + (ir - x_od)**2 < oddistance**2
    roi[circle1] = 1
    roi[circle2] = 0
    plt.imshow(roi,cmap='gray');plt.axis('off');plt.show()

    roi_segment = roi*segment
    #plt.imshow(roi_segment,cmap='gray');plt.axis('off');plt.show()  #显示灰度图像

    centerline,radius=morphology.medial_axis(roi_segment,return_distance=True)
    #plt.imshow(centerline,cmap='gray');plt.axis('off');plt.show()
    #plt.imshow(radius,cmap='gray');plt.axis('off');plt.show()
    vesselcaliber=centerline * radius
    print(vesselcaliber.max())

    #提取特征
    x_temp,y_temp=np.nonzero(vesselcaliber)
    featureset=np.zeros((len(x_temp),3),dtype='float32')
    loc = np.zeros((len(x_temp), 2))
    LABLE=np.zeros((len(x_temp),1))
    t=0
    for i in range(len(x_temp)):
        if vesselcaliber[x_temp[i], y_temp[i]]>0 and vesselcaliber[x_temp[i], y_temp[i]]<10000:
            # 确定标签
            if (label_image[x_temp[i], y_temp[i], 0] > 0):
                TEMP=1
            elif (label_image[x_temp[i], y_temp[i], 2] > 0):
                TEMP=0
            else:
                TEMP=findlabel(label_image, x_temp[i], y_temp[i])
            if TEMP==0 or TEMP==1:
                LABLE[t]=TEMP
                #特征提取
                featureforone = feature_extract1(img,x_temp[i], y_temp[i])
                featureset[t]=featureforone
                loc_temp=np.array([x_temp[i], y_temp[i]])
                loc[t]=loc_temp
                t=t+1

    featureset=featureset[0:t]
    featureset = featureset.astype(float)
    featureset=preprocessing.scale(featureset)#特征归一化
    location=loc[0:t]
    LABLE=LABLE[0:t]
    LABLE=LABLE.ravel()#标签
    print(LABLE.shape)
    x_data.append(featureset)
    y_data.append(LABLE)
    LOCATION.append(location)

np.save('x_data.npy', x_data)
np.save('y_data.npy', y_data)
np.save('LOCATION.npy', LOCATION)





